#pragma once
#include <vector>
#include <array>
#include <memory>
#include <functional>
#include <fstream>
#include <iostream>
#include "SurfaceExtractor.h"
#include "KDTree.h"
#include <Eigen/Dense>
#include <variant>

/**
 * @brief 隐式曲面类
 */
class ImplicitSurface {
public:
    using ImplicitFunction = std::function<double(const Point3D&)>;
    using BoundingBox = std::pair<Point3D, Point3D>;

    ImplicitSurface(ImplicitFunction phi, ImplicitFunction psi = nullptr, 
                   BoundingBox bounds = {{0,0,0}, {0,0,0}})
        : phi_(phi), psi_(psi), bounds_(bounds) {}

    double evaluatePhi(const Point3D& p) const { return phi_(p); }
    double evaluatePsi(const Point3D& p) const { return psi_ ? psi_(p) : 0.0; }
    const BoundingBox& getBounds() const { return bounds_; }

private:
    ImplicitFunction phi_;  // 隐式函数
    ImplicitFunction psi_;  // 边界函数（仅断层使用）
    BoundingBox bounds_;    // 包围盒
};

/**
 * @brief 隐式函数构建器类
 */
class ImplicitBuilder {
public:
    explicit ImplicitBuilder(double gridSize)
        : gridSize_(gridSize) {}

    using FaultSurfaceResult = std::variant<
        std::shared_ptr<ImplicitSurface>,
        std::vector<std::shared_ptr<ImplicitSurface>>
    >;

    /**
     * @brief 构建断层界面隐式函数
     */
    FaultSurfaceResult buildFaultSurface(const std::vector<FaultNode>& faultNodes);

    /**
     * @brief 构建地层界面隐式函数
     */
    std::shared_ptr<ImplicitSurface> buildStrataSurface(
        const std::vector<Point3D>& points);

private:
    /**
     * @brief 加密断层柱线
     */
    std::vector<Point3D> densifyFaultPillars(
        const std::vector<FaultNode>& faultNodes);

    /**
     * @brief 检测是否为特殊平面
     */
    bool detectPlane(const std::vector<Point3D>& points, 
                    char& planeType, double& value);

    /**
     * @brief 检测断层是否为平面
     */
    bool detectFaultPlane(const std::vector<Point3D>& points,
                         Point3D& normal, double& d);

    /**
     * @brief 验证隐式函数
     */
    bool validateSurface(const ImplicitSurface& surface,
                        const std::vector<Point3D>& points,
                        bool isFault = false,
                        const std::string& faultName = "");

    /**
     * @brief 计算点集的法向量
     */
    std::vector<Eigen::Vector3d> computeNormals(
        const std::vector<Point3D>& points,
        bool ensureUpward = true) const;

    /**
     * @brief 构建RBF插值器
     */
    ImplicitSurface::ImplicitFunction buildRBFInterpolator(
        const std::vector<Point3D>& points,
        const std::vector<Eigen::Vector3d>& normals);

    /**
     * @brief 计算点到平面的有向距离
     */
    double computeSignedDistance(
        const Point3D& point,
        const Point3D& planePoint,
        const Eigen::Vector3d& normal) const;

    /**
     * @brief 构建KD树
     */
    std::shared_ptr<KDTree> buildKDTree(const std::vector<Point3D>& points) const;

    /**
     * @brief 构建断层边界函数
     */
    std::function<double(const Point3D&)> buildFaultBoundaryFunction(
        const std::vector<Point3D>& points);

    /**
     * @brief 计算RBF基函数值
     */
    double rbfKernel(double r, double epsilon = 1.0) const;

    /**
     * @brief 计算点到线段的距离
     */
    double pointToSegmentDistance(
        const Eigen::Vector2d& point,
        const Eigen::Vector2d& segStart,
        const Eigen::Vector2d& segEnd) const;

    /**
     * @brief 判断点是否在多边形内
     */
    bool isPointInPolygon(
        const Eigen::Vector2d& point,
        const std::vector<Eigen::Vector2d>& polygon) const;

    /**
     * @brief 添加findKNearestNeighbors的声明
     */
    void findKNearestNeighbors(
        const Eigen::MatrixXd& points,
        int k,
        Eigen::MatrixXi& indices,
        Eigen::MatrixXd& distances) const;

    /**
     * @brief 添加新的辅助方法
     */
    bool isSparseCase(const std::vector<FaultNode>& faultNodes) const;
    std::vector<std::shared_ptr<ImplicitSurface>> 
    buildSparseFaultSurface(const std::vector<FaultNode>& faultNodes);
    std::vector<std::vector<FaultNode>> groupNodesByPillar(
        const std::vector<FaultNode>& faultNodes) const;
    std::vector<std::pair<std::vector<Point3D>, std::vector<Point3D>>> 
    generatePlanePairs(const std::vector<std::vector<FaultNode>>& pillars) const;

    double gridSize_;

    // 完全禁用日志输出
    void debug_log(const std::string&) const {}
    void debug_log(const std::string&, double) const {}
    void debug_log(const std::string&, int) const {}
    void debug_log(const std::string&, size_t) const {}

    // 添加内存预分配的辅助函数
    template<typename T>
    void reserve_vector(std::vector<T>& vec, size_t size) {
        try {
            vec.reserve(size);
        } catch (const std::bad_alloc&) {
            // 静默失败，继续使用默认分配
        }
    }
}; 